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Abstract 

We propose a kinetic Ising model to study phase separation driven 
by surface diffusion. This model is referred to as Model S, and consists 
of the usual Kawasaki spin-exchange kinetics (Model B) in conjunction 
with a kinetic constraint. We use novel multi-spin coding techniques 
to develop fast algorithms for Monte Carlo simulations of Models B 
and S. We use these algorithms to study the late stages of pattern 
dynamics in these systems. 



1 Introduction 



Consider a homogeneous binary (AB) mixture, which is rendered thermody- 
namically unstable by a rapid temperature quench below the miscibility gap. 
The system prefers to be in a phase-separated state at the lower temper- 
ature. The far-from-equilibrium evolution of the system from the unstable 
homogeneous state to the segregated state has received considerable atten- 
tion PI 12 El 13 • This evolution is characterized by the emergence and growth 
of domains enriched in the components A and B. The terms used to describe 
this nonequilibrium process are phase ordering dynamics, domain growth or 
coarsening. A quantitative characterization of phase ordering systems focuses 
on (a) the domain growth law; (b) the statistical properties of the evolution 
morphology; and (c) the temporal correlation of pattern dynamics. 

The equilibrium phase-separated state is uniquely determined by its ther- 
modynamic properties. However, there is a diverse range of kinetic pathways 
which enable segregation. For example, phase separation in alloys is usu- 
ally driven by vacancy-mediated diffusion On the other hand, for fluid 
mixtures, hydrodynamic velocity fields enable convective transport of ma- 
terial along domain boundaries and give rise to novel asymptotic behaviors 
[H]. Furthermore, phase separation in mixtures can be frozen (or near-frozen) 
into mesoscopic states by the presence of quenched disorder [7j, viscoelastic 
effects jHlEl, etc. 

In this paper, we present results from a comparative Monte Carlo (MC) 
study of two kinetic Ising models for phase separation in binary mixtures. 
The first of these is the usual Kawasaki spin-exchange model OE], which 
mimics segregation via diffusion. The second model mimics the case where 
only surface diffusion is permitted. An important goal of this paper is 
methodological, viz., the formulation of a kinetic Ising model where bulk 
diffusion is suppressed. Another important goal is to compare pattern dy- 
namics for phase separation with and without bulk diffusion. 

This paper is organized as follows. In Sec. 2, we describe the kinetic Ising 
models studied here and our MC simulation techniques. Our MC approach 
uses novel multi-spin coding techniques, which enable large-scale and long- 
time simulations of these models. In Sec. 3, we discuss the domain growth 
laws which arise from bulk and surface diffusion, and also present detailed 
numerical results. Finally, Sec. 4 concludes this paper with a summary and 
discussion of our results. 
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2 Numerical Methodology 



2.1 Kinetic Ising Models 



The standard model for binary mixtures is the Ising modeh 



n 



(ji = ±1, 



(1) 



where the spins {uj} {i = 1 ^ N) are located on a discrete lattice. The states 
(Tj = +1 or — 1 denote the presence of an A-atom or B-atom at site i, respec- 
tively. We consider the case with ferromagnetic (J > 0) nearest-neighbor 
interactions, denoted by the subscript (ij) in Eq. (jT)). The phase diagram 
for the binary mixture is obtained in an ensemble with fixed temperature T 
and magnetization M = 

The Ising system does not have an intrinsic dynamics as the Poisson 
brackets (or commutators) of spin variables are identically zero. Therefore, 
one introduces stochastic kinetics by placing the system in contact with a 
heat-bath which induces fluctuations. The Ising model, in conjunction with 
a physically appropriate spin kinetics, is referred to as a kinetic Ising model 
|31 nni- An important example is the Kawasaki spin-exchange model 0, 
which has nearest-neighbor spin exchanges with Metropolis acceptance prob- 
abilities. In an MC simulation of this model, a pair of nearest-neighbor sites 
i and j is randomly selected, and the spins cTj and aj are exchanged. The 
probability that this exchange is accepted is given by 



Here, ATi. is the energy change due to the proposed spin exchange, and 
P = {kBT)~^ is the inverse temperature, with ks denoting the Boltzmann 
constant. In Eq. (0), Li denotes the nearest-neighbors of i on the lattice. A 
single Monte Carlo step (MCS) corresponds to such attempted exchanges. 
A large number of MC simulations of the Kawasaki model have been reported 
in the literature [HI [121 • 

The phase-separation kinetics in this microscopic model is analogous to 
that for the coarse-grained Cahn-Hilliard-Cook (CHC) equation, which is 



P 



min [l,exp(-/3A7^)] , 
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obtained as follows: 



-V ■ J(f, t) 



(3) 



Here, ip{f,t) is the order parameter at space point r and time t. Typically, 
ip{r,t) = pAif.t) — pB{r,t), where pa and ps denote the local densities of 
species A and B. In Eq. the quantities J,D and p denote the current, 
diffusion coefficient, and chemical-potential difference between A and B, re- 
spectively. The chemical potential is obtained as a functional derivative of 
the Helmholtz free energy, which is often taken to have the ip'^-form: 



n-TS 

df 



(4) 



where we have identified (cxj) = '4){fi) in Eq. (^J) and Taylor-expanded vari- 
ous terms. Here, T^. denotes the critical temperature. Finally, the Gaussian 
white noise term 9{f, t) in Eq. (j21) has zero average and obeys the appropriate 
fluctuation-dissipation relation. The CHC equation is also known as Model 
B in the Hohenberg-Halperin classification scheme for critical dynamics • 
Further, using a master-equation approach, the CHC equation can be mo- 
tivated from the spin-exchange model |l4j. Therefore, we will subsequently 
refer to the Kawasaki model as "Model B". 

Before proceeding, we stress that the general form of the CHC equation 
contains an order-parameter-dependent mobility I16| IT7j: 



D{^) = Do f 1 - ^ 



(5) 



where Vo is the saturation value of the order parameter at T = 0. This 
is not consequential for quenches to moderate temperatures, but plays an 
important role for deep quenches where ip ~ ±-00 in bulk domains. In that 
case, bulk diffusion is effectively eliminated and domain growth proceeds by 
surface diffusion [T8 t ll9 | l2nj. In the context of the Kawasaki model, this can 
be understood by focusing on an interfacial pair with the minimum barrier 
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for interchange: (Jj = +1 at the periphery of an up-rich domain and having 
only one neighbor with the same spin value, and o"j = — 1 in a down-rich 
domain. At low temperatures, the bulk domains are very pure and the energy 
barrier to the interchange 0", aj is A7i = 4J. Thus, the time-scale for 
this interchange tk ~ exp(/3A7i) ^ oo as T ^ 0, effectively blocking bulk 
diffusion. Of course, once an impurity spin is placed inside a bulk domain, 
there is no further barrier to its diffusion. 

Apart from this natural blocking of bulk diffusion at T = 0, there are 
systems where the bulk mobility diminishes drastically due to physical pro- 
cesses, e.g., one or both of the components may undergo a glass or 
gelation [221 12^] transition. At the phenomenological level, this has been 
modeled by setting the mobility to zero in regions rich in the glass-phase or 
gel-phase. At the microscopic level, we propose a kinetic Ising model where 
bulk diffusion is suppressed by introducing a kinetic constraint. We disallow 
exchanges ^ aj if the neighboring spins of the pair are all parallel, even 
though such an exchange would not raise the energy. In this case, segrega- 
tion is driven primarily by diffusion along domain boundaries, though some 
bulk transport occurs via impurity n-spin clusters. (This bulk diffusion is 
negligible for moderate to deep quenches.) We will subsequently refer to this 
model as "Model S" [TH]. Clearly, Model S can be generalized to the case 
of reduced (though non-zero) mobility in the bulk domains. This is done by 
allowing spin-exchanges with different time-scales depending on the number 
and type of parallel neighbors for a spin pair. 

2.2 Numerical Details 

All our MC simulations were performed on an L x L square lattice with 
periodic boundary conditions. At time t = 0, the temperature was quenched 
from T = oo to T < Tc, where Tc ~ 2.269 is the critical temperature of the 
d = 2 Ising model. (All energy scales are measured in units of J, and we set 
the Boltzmann constant ks = 1.) The disordered initial state consisted of a 
uniform mixture of A^^ A-atoms and Nb B-atoms with = A^^ + Nb- The 
case with A^^ = Nb corresponds to a critical quench. 

Our MC simulations exploit the technique of multi-spin coding. For a 
general introduction to this technique, see Ref. [21]. The basic idea is that 
one can exploit the 64-bit computer architecture to undertake a parallel sim- 
ulation of 64 systems. This is done by storing the spin cTj at site i in the k^^ 
system in the k^^ bit of a 64-bit word S[i]. Recall that, in one elementary 
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move for Model B, we propose to exchange the spins located on nearest- 
neighbor sites i and j. For Model S, we impose the kinetic constraint that a 
pair of spins surrounded by aligned neighbors is not exchanged. 

For the d = 2 square lattice considered here, each site has four nearest- 
neighbors. Let no, rii and n2 be the three nearest-neighbors (other than 
j) of site i. Similarly, let mo, mi and m2 be the three nearest-neighbors 
(other than i) of site j. To determine the change in energy resulting from 
the proposed spin exchanges in all 64 simulations, we first identify which of 
the six neighbors (no, ni, n2, mo, mi, m2) are antiparallel. This can be done 
in six operations with the exclusive or operation ©: 

Ak = S[i]®S[nk], k = 0^2, 

Bk = S[j]®S[m,l k = 0^2. (6) 

The energy change associated with the spin exchange and (thereby) the ac- 
ceptance probability is governed by the number of antiparallel spins. In 
an ordinary program, this would involve a summation over the surrounding 
spins. With logical operations, it is more convenient to determine the logical 
variables Pk which tell whether o"j is antiparallel to at least k of its neighbors 
(other than (Xj). Note that the Metropolis algorithm only requires Pi, P2 and 
P3. These can be obtained with six operations: 



P2 


= AoAAi, 


Pi 


= AoVAi, 


P3 


= A2AP2, 


P2 


= P2V(A2APi), 


Pi 


= P1VA2. 



Similarly, the variables Qk that tell whether aj is antiparallel to at least k of 
its neighbors (other than cTj) can be obtained with six operations. 

Finally, the acceptance probabihty for the proposed spin exchanges is 
obtained by using random bit patterns Rq, Ri and i?2- These are designed 
so that the probability for each bit to be 1 is Pf, = exp(— 4/3J). Thus, the 
following statements comprise the core of our Model S algorithm: 

Flip = (5[^]©5[j]) A(Pi vgsVPo) A(P2Vg2VPi)A(P3VQiVP2), 

Flip = FlipA(Pi V-Q3) A(Qi V-P3), 

S\i] = 5[2]©Flip, 

S[j] = ^b]©FliP- (8) 
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The implementation of Model B dynamics is simply obtained by omitting 
the second of the above statements. 

These 36 operations for Model S (or 30 for Model B) act on all 64 bits 
and thus perform 64 elementary moves. In conjunction with the required 
load and store operations, and generation of the random bit patterns, our 
implementation of Model B for a 512^ system with multi-spin coding requires 
4.6 ns CPU-time per elementary move on an AMD-64 computer with 3 GHz 
clock frequency. This should be contrasted with a direct implementation of 
this model, which requires approximately 100 ns CPU-time per elementary 
move on the same machine. 

The procedure outlined above, which simulates 64 separate systems, is 
known as a synchronous multi-spin algorithm [21] . The boundaries of these 
64 systems can be glued together to yield an asynchronous multi-spin algo- 
rithm simulating one system which is 64 times larger. This comes at the 
cost of (a) more complicated programming; and (b) a small reduction in the 
program efficiency. The statistical results for domain morphologies presented 
in Sees. 3.1, 3.2 and 3.3 were obtained by averaging over 150 asynchronous 
simulations with system sizes L = 512. The results for the autocorrelation 
function in Sec. 3.4 were obtained by averaging over 64 synchronously simu- 
lated systems with L = 1024. 

3 Detailed Results 

As stated earlier, the initial condition for our MC simulations consists of a 
random configuration. The temperature is quenched to T < Tc at t = 0, and 
the system evolves via either Model B or Model S dynamics towards its new 
equilibrium state. Figure Q shows the typical time evolution for a critical 
composition (50 % A and 50 % B) after a quench to T = 0.63Tc for Model 
B (left) and Model S (right). Notice that the evolution morphology has a 
characteristic domain size, which we denote as R{t). The growth process 
is substantially slower for S-dynamics, as expected. We will demonstrate 
shortly that the growth law due to bulk diffusion is R{t) ~ t^^^, which is 
referred to as the Lifshitz-Slyozov (LS) growth law. The corresponding law 
for segregation via surface diffusion is R{t) ~ t^^^. However, we reiterate that 
bulk diffusion is not eliminated entirely in Model S because of the presence 
of impurity spin clusters in bulk domains. At high temperatures, there is a 
reasonable fraction of impurity spins and we expect the S-dynamics to cross 
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Figure 1: Evolution pictures for phase separation in a binary (AB) mixture 
with a critical composition. The component A (cr = +1) is marked in black, 
and the component B (cr = —1) is unmarked. The system was quenched at 
time t = from T = oo to T = 0.63Tc. The top and bottom panels show 
snapshots at times t = 10^ and 10^ MCS, using either Model B dynamics 
(left), or Model S dynamics (right). The MC simulations were done on square 
lattices of size 512^ with periodic boundary conditions. The details of the 
simulations are described in the text. 

over to t^/^-growth at late times. The crossover time increases rapidly at 
lower temperatures where there are very few impurity spins in the bulk as 
Pimp ^ [1 + exp(8/?J)]"^ md = 2. 

We will study the evolution depicted in Fig. [T] using quantities like the 
correlation function and autocorrelation function. 

3.1 Growth Laws 

The first relevant property is the growth law governing the segregation pro- 
cess. We computed the typical domain size R{t) as the first zero-crossing of 
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the two-point correlation function: 



N 



C{r,t) 



1=1 



= 9 



R 



(9) 
(10) 



Here, r denotes the displacement vector, and we consider systems which are 
translationally invariant and isotropic. The angular brackets in Eq. ^ de- 
note an averaging over independent initial conditions and noise realizations. 
Equation (fTUj) is the dynamical-scaling property of the correlation function 
and reflects the fact that the morphology is self-similar in time, upto 
a scale factor (see Fig. Q). One can use other definitions of the length scale 
also, but these are all equivalent in the scaling regime. 

At this stage, it is useful to clarify the domain growth laws which arise 
due to bulk and surface diffusion. A convenient starting point is the CHC 
equation Q with an order-parameter-dependent mobility D^ip). We consider 
a general situation where the diffusion coefficient at the interface [ip = 0) is 
Dg, and that in the bulk [i/j = V's(^)] is -Dt with Dh < Dg. This difference 
in surface and bulk mobilities may be the consequence of low-temperature 
dynamics or due to physical processes like glass formation or gelation. Then, 
the simplest functional form which models the mobility is 



which is equivalent to Eq. (0) with Do 
the deterministic case of Eq. Q: 



a 



(11) 



Ds and ipQ = ipl/a. We focus on 



^^(r,t) 



D.^ ■ 



a- 



^2 



-(r,-T)^ + ^^^^ 



(12) 



where we have used the ip'^-foim of the free energy from Eq. Q. The sat 
uration value of the order parameter in Eq. (|12j) is ips{T) = 



3(1-T/T,). 

Using the natural scales for the order parameter, length and time, we can 
rewrite Eq. (fT^ in the dimensionless form: 



d_ 
di' 

a G [0, 1] for D^, < D^ 



^(f,t) = V- ( 1 - V (-^/' + ^/-^ - VV 



(13) 
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The RHS of Eq. (fT5|) can be decomposed as [201 



[-^ + ^^ -V^ip)], (14) 



where the first term on the RHS corresponds to bulk diffusion. This term 
disappears for a = 1 or = 0. The second term on the RHS corresponds 
to surface diffusion and is only operational at interfaces where ~ 0. Fol- 
lowing Ohta |2Sl! we can obtain an equation which describes the interfacial 
dynamics. The location of the interfaces fi{t) is defined by the zeros of the 
order-parameter field: 

^[rl(t),t] =0. (15) 

Focus on a particular interface, and designate the normal coordinate as n 
(with dimensionality 1) and the interfacial coordinates as a [with dimension- 
ality (d—l)]. Then, the normal velocity Vn{a, t) obeys the integro-differential 
equation [^1^: 



4 I da'G[fi{a),fi{a')]vn{a',t) ~ {1 - a)aK{d,t) + 

4a / da'G[fi{d),fi{a')]V'^K {a', t), (16) 



where K{d, t) is the local curvature at point a on the interface, and a is the 
surface tension. The Green's function G(x, x') obeys 

- V^G(f,x') = 5(x-x')- (17) 

A dimensional analysis of Eq. (jlfij) in the scaling regime yields the growth 
laws due to surface and bulk diffusion. We identify the scales of various 
quantities in Eq. (fT^ as 

[da] ~ R'^-\ [G] ~ R^~\ 

, [K]^R-\ (18) 

This yields the crossover behavior of the length scale as 

R{t) ~ (a^)^/^ t<tc, 



[(1 - a)atY'\ t > t„ (19) 
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where the crossover time is 



a 

1 -«)V4' 



(20) 



The above scenario apphes for both Models B and S, as < in either 
case. At moderate temperatures, this crossover occurs rapidly for Model B in 
our simulations. However, in Model S, there is a drastic suppression of bulk 
diffusion with Di, <^ Dg and a ~ 1. Therefore, the crossover to t^/^-growth 
is strongly delayed and not observed over simulation time-scales. This is 
seen in Fig. which plots R vs. t at temperatures T = 0.63Tc and O.SSTc 




10 
t [MCS] 



10° 10" 



10=^ 
t [MCS] 



Figure 2: Typical domain size [R) as a function of time (t) after a quench 
at t = from T = 00 to T = 0.63T, (left) and 0.88T, (right). The circles 
and squares indicate data obtained with Model B and Model S dynamics, 
respectively. Lines with slope 1/4 and 1/3 are also provided on the plots as 
guides to the eye. 



for Models B and S. The data for Model S is consistent with the growth law 
R ~ t^/^. The early-time data for Model B is also consistent with this growth 
law, as surface diffusion is dominant at early times. At late times, one sees 
crossover behavior between the t^/'^-regime and the asymptotic t^/^-regime. 
Note that the crossover for Model B is delayed at the higher temperature 
T = 0.88Tc because the decrease in a is more than compensated by the 
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reduction in a due to softening of the interfaces as T — > T~ [see Eq. ((201)] • 
More generally, we stress that it has been notoriously difficult to observe the 
asymptotic t^/^-growth in MC simulations of the Kawasaki model fTJ ^] . 
Similar results have been obtained from Langevin studies of coarse-grained 
models [IHl CHI 1201 ■ 

Before proceeding, we remark that we have also studied models where 
bulk diffusion is more strictly suppressed by imposing additional kinetic con- 
straints which eliminate 2-spin diffusion, 3-spin diffusion, etc. The corre- 
sponding results for the growth law are numerically indistinguishable from 
the Model S results in Fig. |21 over the time-scales of our simulation. This un- 
derlines the utility of the proposed Model S in the context of phase separation 
via surface diffusion. 

It is also relevant to discuss off-critical quenches, where one of the com- 
ponents is present in a larger fraction. In Fig. El we show evolution pictures 
for Models B and S for the case with 25% A and 75% B. If the evolution 



Figure 3: Analogous to Fig. ^ but for an off-critical quench with 25% A and 
75% B. The temperature is T = 0.63Tc. 
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morphology is not bicontinuous, e.g., there are droplets of the minority phase 
in a matrix of the majority phase, the surface- diffusion mechanism is unable 
to drive growth. Nevertheless, at high temperatures, growth may still pro- 
ceed by the Brownian motion of droplets . The corresponding growth law 
depends explicitly on the dimensionality: 

R{t) - (Tt)i/('^+2)_ ^21) 

Thus, domain growth through droplet motion obeys the law R{t) ~ (TtY^^ 
in d = 2, which is analogous to the surface-diffusion growth law. The growth 
kinetics of Models B and S for off-critical mixtures at T = 0.63Tc is shown in 
Fig. I3J We see that the S-dynamics shows the expected t^^^-growth over ex- 
tended time-regimes. As before, the B-dynamics shows a crossover behavior 
between the t^/^-regime and the asymptotic t^/^-regime. At low tempera- 
tures, the Brownian mechanism is ineffective and the evolution of Model S 
freezes into a meso-structure. 

3.2 Correlation Functions 

Next, let us study the morphological features of the evolution in Figs. [T] and 
El These are usually characterized by (a) the correlation function defined in 
Eq. 0, or (b) its Fourier transform, the structure factor. We have confirmed 
that these quantities exhibit dynamical scaling for both Models B and S. For 
the sake of brevity, we do not show these results here. 

An important theme in this context is a comparison of the morphologies 
arising from both dynamics. Earlier studies with coarse-grained models JHl 

1201 have found that the correlation functions and structure factors are 
numerically indistinguishable for growth driven by bulk diffusion or surface 
diffusion. At the visual level, this also seems to be suggested by the snapshots 
in Figs.^andEl In Fig.^l we compare the scaling functions for Models B and 
S for a critical quench with T = 0.63Tc. To eliminate finite-size effects, we 
consider cases with the same typical domain size: t = 10^ MCS in Model S 
and t = 3.4 X 10^ MCS in Model B. In Fig.E^a), we plot C(r, t) vs. r/R. The 
scaling functions superpose on the scale of the plot, in accordance with earlier 
studies of phenomenological models. In Fig. El^b), we plot C{r,t) ■ {r/R) vs. 
r/Rso that the large-distance behavior is magnified. Some subtle differences 
between the two functions are seen at large distances r/R > 2. We make 
some observations in this regard: 
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Figure 4: Analogous to Fig. |2l but for an off-critical quench with 25% A and 
75% B. We show results for T = 0.63T^. 

(a) The statistical error in the difference between the curves at the second 
peak is about four times smaller than the difference, so it cannot be attributed 
to statistical fluctuations. 

(b) We have also replotted the correlation functions for Models B and S 
from different times on the scale C{r,t) ■ {r/R) vs. r/R. In that case, the 
secondary peaks show a much better collapse, suggesting that the discrepancy 
in Fig. E^b) is not the result of corrections to scaling. 

Though it is difficult to attribute physical significance to these differences, 
it is conceptually important to stress the observable differences between the 
morphologies for Models B and S. Similar scaling plots for the off-critical 
quench shown in Fig. El are shown in Fig. El It is known that the scaling 
functions for phase-separating systems depend on the degree of off-criticality 



14 



012345 012345 
r/R r/R 



Figure 5: Superposition of scaling functions for Models B (solid line) and 
S (dashed line) for the evolution depicted in Fig. [T] For Model S, the data 
set corresponds to t = 10^ MCS; for Model B, the data set corresponds to 
t = 3.4 X 10^ MCS. Both domain sizes coincide at these times, (a) Plot of 
C{r,t) vs. r/R. (b) Plot of C{r,t) ■ {r/R) vs. r/R, so as to magnify the tail 
behavior. 

|28j . Notice that the oscillations in the plot of C{r,t) vs. r/R diminish with 
increase in the off-criticality. Further, the discrepancy between the scaling 
functions for Models B and S is larger for the off-critical case. 

3.3 Island Distribution and Excess Energy 

Our subsequent results will focus on the case of a critical quench. An alter- 
native method of describing the domain morphology is the island-size distri- 
bution. We define an island as a set of aligned spins, all of whose neighbors 
are either part of the island, or have an antiparallel spin. Tafa et al. 2^ 
have shown that the domain-size distribution in a phase-separating system 
P{1, t) exhibits scaling, and has an exponential decay: 



where / is the domain size and a is a constant. The corresponding scaling 
form for the island-size distribution P{s,t) in d = 2 is obtained as 




for X 



oo. 



(22) 
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Figure 6: Analogous to Fig. but for an off-critical quench with 25% A and 
75% B at T = 0.63Tc. In this case, the times for the different data sets are 
t = 9.8 X 10^ MCS (Model S) and t = 2.7 x 10^ MCS (Model B). 



= 3{x) = ^f{V^), (23) 

where 6 is a geometric factor, and (s) is the average island size. 

In Fig. HJ we plot ,/s{s)P{s, t) vs. ^s/{s) for both Models B and S. We 
make two observations in this context. First, the data for the two models is 
numerically indistinguishable on the scale of this plot. The subtle differences 
in the correlation-function data are not seen in the island-size distribution 
function. Second, the plot in Fig.[7|exhibits an exponential decay, as expected 
from Eqs. 

A macroscopic quantity which depends on the density of small islands 
is the total energy E{t). The interfacial energy for a domain is aW^'^, and 
the number of domains in the system ~ N/R^. Thus, the overall interfacial 
energy depends on the length scale as E{t) — E{oo) ~ Na/R. In Fig. |Hl we 
plot E{t)/N vs. R-^ for both Models B and S at T = 0.63Tc and O.SST,. We 
observe a power-law convergence of the excess energy with the slope being 
proportional to the surface tension cr(T). Again, the data sets at the same 
temperature cannot be distinguished on the scale of the plot. 
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Figure 7: Scaled probability distributions for island-sizes in Models B (solid 
line) and S (dashed line) at T = 0.88Tc. The data is shown on a linear-log plot 
as ,/s{^P{s,t) vs. ^s/{s), suggested by Eqs. ((12))- ((221) • For Model S, the 
data set corresponds to t = 10® MCS; for Model B, the data set corresponds 
to t = 3.4 X 10^ MCS. Both domain sizes coincide at these times. 



3.4 Aging of the Autocorrelation Function 

The data presented so far has focused on the morphological features of the 
phase-separating system. Let us next study the temporal correlation of the 
pattern dynamics in Fig.^ This is measured by the autocorrelation function: 

1 ^ 

A{t^,T) = —J2[{a^it^)aiit^ + T))-{ai{t^)){ai{t^ + T))], (24) 
i=i 

where the times and (t^ + r) are measured after the quench at t = 0. Here, 
is the reference time for measurement of the autocorrelation function, and 
is referred to as the waiting time. The most general correlation function cor- 
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Figure 8: Plot of the energy per site E{t)/N vs. for Model B at tem- 
peratures T = 0.63Tc (squares) and T = 0.88Tc (circles), and Model S at the 
same two temperatures (diamonds and crosses). 

responds to unequal space and time, and combines the definitions in Eqs. Q 
and (|24p . Equilibrium systems are stationary and the corresponding A{tw, r) 
only depends upon the time-difference r. On the other hand, for nonequilib- 
rium systems, A{tyj, r) depends on both t^j and r. 

There have been some earlier studies of A{tw, t) for domain growth in ki- 
netic Ising models. There are two mechanisms which drive the decorrelation 
process: 

(a) First, there are fluctuations in bulk domains, which give a stationary 
contribution. Huse and Fisher jSQ, studied decorrelation arising from the ap- 
pearance of a droplet of (say) down-spins in an up-domain. The probability 
that a droplet of size R appears via fluctuations is Pd oc exp(— /3cri?'^~^). The 
lifetime of this droplet is r ~ i?^/*^, where is the growth exponent. Thus, 
the corresponding autocorrelation function shows a stretched-exponential be- 
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havior: 



Aq(r) ~exp(-/5ar^), 

6 = {d — 1)0 for d < dc and 1 for d > d^. 



(25) 



Here, the critical dimensionality is defined by {dc — 1)0 = 1- 
(b) Second, there is decorrelation due to domain-wall motion. This can 
be either stochastic (due to thermal fluctuations) or systematic (due to the 
curvature- reduction mechanism). Consider the T = case, where there are 
no fluctuations in the bulk or the surface. The characteristic domain-wall 
velocity decreases with time, so this mechanism gives a non-stationary or ag- 
ing (tu;-dependent) contribution [21]. Fisher and Huse [S2] used scaling ideas 
to argue that the aging contribution to A(t^, r) has a power- law dependence 
on the length scale: 



There have been various studies of the aging exponent A in cases with 
both spin- flip and spin-exchange kinetics . For power-law domain growth, 
Eq. obeys the scaling form A^g^i'tw^'T) = h^r/tu,), which has been ob- 
served in some studies of spin glasses pi] . 

In Fig. El we plot A(tw, t) vs. r for Models B and S for a critical quench 
to T = 0.63Tc. The solid lines denote data for Model B with waiting times 
tu, = 10^, 10^, 10^, 10^ (from left to right). The dashed lines denote the corre- 
sponding data for Model S. As expected, the autocorrelation function decays 
more rapidly for Model B. We make the following observations in this regard: 

(a) The quantity y4(tu,,r) exhibits aging, with an explicit dependence on t^, 
for both Models B and S. In general, the decay is slower for larger t^, i.e., 
when the domain size of the reference state is larger. Further, the decay is 
faster for higher temperatures, where larger fluctuations are present in the 
system. 

(b) The data in Fig. El is plotted on a log-log scale, and exhibits a contin- 
uous curvature for both Models B and S. This is not consistent with the 
simple power-law decay in Eq. (f^ . As a matter of fact, the autocorrelation 
data does not even exhibit r/t^-scaling, as we have confirmed. In the case 
of Model B, this is because the decorrelation process is driven by both bulk 
fluctuations (with a stationary contribution) and domain-wall motion (with a 
non- stationary contribution). In the case of Model S, bulk fluctuations have 




(26) 



R{t^ + t) 
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Figure 9: Time-dependence of the autocorrelation function for Models B 
(solid line) and S (dashed line) at T = 0.63Tc. The waiting times are = 
10^, 10^, 10^ 10^ MCS (from left to right). 

been effectively eliminated and one may naively expect to recover power-law 
decay. However, this is not the case because interfacial fluctuations also con- 
tribute to decorrelation. We believe that the scaling behavior in Eq. ()26p is 
only realized in kinetic Ising models or their coarse-grained analogs at T = 0. 
In this limit, coarsening occurs only through the systematic motion of inter- 
faces and the system always reduces its energy. However, the T = limit 
is not interesting in the context of kinetic Ising models because the evolving 
system invariably gets trapped in local free-energy minima, 
(c) In recent work, Puri and Kumar have studied the decorrelation pro- 
cess in a spin-1 model using a stochastic model based on the continuous-time 
random walk formalism. We are currently trying to adapt their modeling to 
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understand the behavior in Fig. |H1 



4 Summary and Discussion 

Let us conclude this paper with a summary and discussion of the results pre- 
sented here. We have studied phase separation in a kinetic Ising model for 
phase separation mediated by surface diffusion. This model (referred to as 
Model S) is obtained by imposing a kinetic constraint on the usual Kawasaki 
kinetic Ising model (referred to as Model B). In general, the surface diffusion 
mechanism can drive segregation only when the morphology consists of per- 
colated clusters, i.e., for near-critical quenches. We have undertaken Monte 
Carlo (MC) simulations of Models B and S using multi-spin coding tech- 
niques. These provide accelerated algorithms which enable the simulation of 
large systems for extended times. Our results show that the major difference 
between the morphologies of Models B and S lies in the growth dynamics. 
In this regard, it is relevant to emphasize the following: 

(a) The early-time dynamics {t <^t^) of Model B is also dominated by sur- 
face diffusion with the growth law R{t) ~ t^/^. For late times (t ^ ), 
there is a crossover to the t^/'^-growth regime. In the limit of T — > 0, the 
crossover time diverges (tf —>■ oo). However, the low-temperature dynamics 
of Model B usually freezes into metastable states. Therefore, it is hard to 
see an extended regime of t^/^-growth in Model B. 

(b) Our kinetic constraint eliminates single-particle bulk diffusion, and we see 
extended regimes of growth driven by surface diffusion. However, n-particle 
diffusion (with n > 2) is still possible and is governed by the probability 
for existence of impurities in bulk domains. Thus, at sufficiently large times 
{t ^ tf), we again expect a crossover to t^/^-growth. However, this crossover 
is extremely delayed, even at moderate temperatures. 

(c) We have also studied kinetic models with constraints which eliminate the 
diffusion of ra-spin clusters. The domain growth data obtained from these 
models is numerically indistinguishable from that for Model S. 

(d) For highly off-critical quenches, the morphology consists of droplets of 
the minority phase in a matrix of the majority phase. In this case, the sur- 
face diffusion mechanism cannot drive phase separation. However, Brownian 
motion and coalescence of droplets also gives rise to t^/^-growth in d = 2. 

Apart from growth laws, we have also studied quantitative properties 
of the evolution morphology like correlation functions and island-size dis- 
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tribution functions. There are subtle differences in the scaled correlation 
functions for Models B and S, but it is difficult to attribute physical signifi- 
cance to these. Further, these differences are not reflected in the island-size 
distribution function. 

Finally, we have studied the aging of the autocorrelation function A(tw, r) 
in Models B and S. In both cases, we find that the decorrelation process is 
driven by both fluctuations and domain-wall motion. Thus, A{t^, r) does 
not exhibit a simple power-law decay or scaling behavior. We are presently 
adapting the continuous-time random walk approach developed by Puri and 
Kumar to study the aging of the autocorrelation functions in Models B 
and S. 
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